#include "libmesh/quadrature.h"
#include "Assembly.h"
#include "FluxJumpIndicator.h"
#include "CFDProblem.h"

template<>
InputParameters validParams<FluxJumpIndicator>()
{
  InputParameters params = validParams<InternalSideIndicator>();
  params += validParams<TransientInterface>();
  params += validParams<MultiVariableInterface>();
  params.addParam<Real>("scale", 1, "人工粘性放大尺度");
//  params.set<VariableName>("variable") = "density";
  return params;
}


FluxJumpIndicator::FluxJumpIndicator(const InputParameters &parameters) :
	InternalSideIndicator(parameters),
    TransientInterface(this),
	MultiVariableInterface(parameters),
	_cfd_problem(static_cast<CFDProblem&>(_fe_problem)),
	_cfd_data(_cfd_problem),
	_cfd_data_neighbor(_cfd_problem),
    _current_elem_volume(_assembly.elemVolume()),
    _neighbor_elem_volume(_assembly.neighborVolume()),
	_current_side_volume(_assembly.sideElemVolume())
{
	// foreach(MooseVariable *var, _variables)
	for(auto var : _variables)
		addMooseVariableDependency(var);
}

void FluxJumpIndicator::computeIndicator()
{
	Real sum = 0;
	for (_qp=0; _qp<_qrule->n_points(); _qp++)
		sum += _JxW[_qp]*_coord[_qp]*computeQpIntegral();

	sum /= _current_side_volume;
	Real scale = getParam<Real>("scale");
	sum *= scale;

	{
		Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
		_solution.add(_field_var.nodalDofIndex(), (sum));
		_solution.add(_field_var.nodalDofIndexNeighbor(), (sum));
	}
}

Real FluxJumpIndicator::computeQpIntegral()
{
	for (size_t i = 0; i < _uh.size(); ++i)
	{
		_cfd_data.uh[i] = (*_uh[i])[_qp];
		_cfd_data_neighbor.uh[i] = (*_uh_neighbor[i])[_qp];
	}
	_cfd_data.reinitInviscous();
	_cfd_data_neighbor.reinitInviscous();

	Real dp_du[] = {
			(0.5*_cfd_data.vel.norm_sq() + 0.5*_cfd_data_neighbor.vel.norm_sq()),
			-(_cfd_data.vel(0)+_cfd_data_neighbor.vel(0)),
			-(_cfd_data.vel(1)+_cfd_data_neighbor.vel(1)),
			-(_cfd_data.vel(2)+_cfd_data_neighbor.vel(2)),
			(1+1)
	};

	Real indicator(0);
	for (int eq = 0; eq < 5; ++eq)
	{
		Real flux_jump= (_cfd_data.invis_flux[eq] - _cfd_data_neighbor.invis_flux[eq]) * _normals[_qp];
		indicator += fabs(flux_jump*(dp_du[eq]));
	}

//	Real p_inf = 1./_cfd_problem._gamma/_cfd_problem._mach/_cfd_problem._mach;
	Real p = fabs(_cfd_data.p + _cfd_data_neighbor.p);// + p_inf;//_cfd_problem._gamma;
	return (indicator)/p*_cfd_problem._gamma;
}

void FluxJumpIndicator::finalize()
{

}

